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m Abstract 
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Independent Component Analysis (ICA) is a computational technique for revealing la- 

X> 

tent factors that underlie sets of measurements or signals. It has become a standard tech- 

PLh 

OO nique in functional neuroimaging. In functional neuroimaging, so called group ICA (gICA) 

,——1 seeks to identify and quantify networks of correlated regions across subjects. This paper re- 

ports on the development of a new group ICA approach, Homotopic Group ICA (H-glC A) , 
^ for blind source separation of resting state functional magnetic resonance imaging (fMRI) 

data. Resting state brain functional homotopy is the similarity of spontaneous fluctuations 

r-H 

^ between bilaterally symmetrically opposing regions (i.e. those symmetric with respect to 

m 

the mid-sagittal plane) (Zuo et al. 20101. The approach we proposed improves network 
estimates by leveraging this known brain functional homotopy. H-gICA increases the po- 
tential for network discovery, effectively by averaging information across hemispheres. It is 



X 



theoretically proven to be identical to standard group ICA when the true sources are both 
perfectly homotopic and noise-free, while simulation studies and data explorations demon- 
strate its benefits in the presence of noise. Moreover, compared to commonly applied group 
ICA algorithms, the structure of the H-glCA input data leads to significant improvement 
in computational efficiency. A simulation study comfirms its effectiveness in homotopic, 
non-homotopic and mixed settings, as well as on the landmark ADIID-200 dataset. From 
a relatively small subset of data, several brain networks were found including: the visual, 
the default mode and auditory networks, as well as others. These were shown to be more 
contiguous and clearly delineated than the corresponding ordinary group ICA. Finally, in 



addition to improving network estimation, H-gICA facilitates the investigation of functional 

homotopy via ICA-based networks. 
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1 Introduction 



The function of the human brain during rest can be investigated using various functional 



measurement techniques (Biswal et aL 1995; Gusnard et aL 2001 Gao et aL, 2009). 
Patterns of resting-state brain functional across individuals provide insights into baseline 
activity of the human brain in the absence of experimental stimuli. 

Resting state functional magnetic resonance imaging (rs-fMRI), obtained using blood 
oxygen level dependent (BOLD) signals, is a key driving force for investigating baseline 
fluctuations in brain activity. Much recent attention has been focused on identifying 
the possible origins and clinical manifestations of variation in the BOLD signals from 
rs-fMRI data. Neuroimaging studies have identified associations between resting state 
brain networks estimated via fMRI data with aging, cognitive function and neurological 



and psychiatric disorders (Damoiseaux et al. 2008; Rombouts et al. , 2005). One popular 
approach for locating putative networks is blind source separation, which decomposes 
neuroimaging data into an outer product of spatial maps multiplied by their respective 
time courses. Notably, blind source separation does not require a specific fMRI paradigm, 
so is applicable to resting state data. Two popular exploratory data analysis techniques 
for blind source separation are Principal Component Analysis (PGA) and Independent 
Component Analysis (IGA). IGA can be distinguished from PGA by its focus on model- 
level independence and non-Gaussianity. Moreover, IGA as a matrix decomposition does 
not yield decomposition vectors that are orthonormal. Finally, IGA is usually applied 
after PGA-based dimensional reduction, and thus can be thought of as a non-orthonormal 
reorganization of PGA. 

IGA has become popular for analyzing neuroimaging data, having been successfully 



applied to single-subject analysis (Guo and Pagnoni 2008 Beckmann and Smith, 2004 



McKeown et al. , 1997). The extension of IGA to group inferences provides common 



independent components across subjects and enables identification of putative underlying 
brain networks for the group. Several multi-subject ICA approaches have been proposed: 



Calhoun et al. (2001b) presented (and named) group ICA (gICA) and created the Matlab 



toolbox (GIFT), which provides estimation. GIFT consists of a first-dimension reduction 
using PCA for each subject, followed by a temporal concatenation of the reduced data. 



after which ICA is then applied to the aggregated data. More recently, Beckmann and 



Smithj ( 2005[ ) proposed a tensor ICA (TIC A) by extending the single-session probabilistic 



ICA (PICA) (Beckmann and Smith, 2004). TICA factors the multi-subject data as a 
tri-linear combination of three outer products, which represents the different signals and 
artifacts present in the data in terms of their temporal, spatial, and subject-dependent 
variations. Other grouping methods proposed perform ICA on each subject and combine 



the output into a group by using self-organizing clustering (Esposito et al. , 2005 ) or spatial 



correlation (Calhoun et al. 2001a). 

Among existing methods, GIFT is perhaps the most commonly used for performing 
group ICA analysis of multi-subject fMRI data. The spatial independence assumed by 
GIFT, and any other spatial ICA algorithms, is well suited to the spatial patterns seen in 



fMRI (McKeown et al. , 1997). Moreover, its empirical performance has been consistently 



vahdated (Sorg et al. 2007 Garrity et al. 2007; Sambataro et al. 2010). 

This paper describes a novel modified group ICA method - homotopic group ICA (H- 
glCA) - to identify the underlying patterns of brain activity for functionally homotopic 
putative networks. Left to right bilaterally symmetric patterns of interhemispheric syn- 
chrony and co-activation are among the most frequent findings in neuroimaging studies 



(Toro et al. , 2008). Functional homotopy, the similarity of spontaneous fiuctuations be- 
tween bilaterally symmetrically opposing regions with respect to the mid-sagittal plane is 



a "fundamental characteristic of the brain's intrinsic functional architecture" ( Zuo et al. 
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2010 ). Via H-gICA, the information of homotopic brain function is utilized to improve the 



identification of underlying brain networks. A spatial independence assumption is made 
relating to all voxels in each hemisphere. In a simulation study, H-gICA is shown to be 
preferable to (our implementation of) GIFT when the actual signals are homotopic and is 
competitive with GIFT under non-homotopic signal conditions. The efficacy of H-gICA 



methodology is demonstrated by an application to the ADIID-200 dataset (Milham, 2012 



Eloyan et al. , 2012). From the fifteen components produced by H-gICA, several common 



brain networks were found, being clearly represented in smoother, more clearly delineated 
and contiguous volumes than ordinary glCA. The main networks found include the visual, 
default mode and auditory. In addition to improving network estimation, H-glCA allows 
for the investigation of functional homotopy via ICA-based networks. Here the quantifi- 
cation of functional homotopy of such networks are defined using a similar concept to that 



of Joel et al. (2011 ). Specifically, by having left and right hemispheric terms in the model, 
investigations of ICA weight matrices includes homotopic information, unlike ordinary 
group ICA. 

The remainder of the paper is organized as follows. ICA and fast, fixed-point algo- 
rithm known as FastICA proposed by (Hyvarinen, 1999) is reviewed in Section |2| and 
subsequently used in the estimation of H-glCA. Section [3] is the theoretical body of the 



paper: Section 3.1 discusses about data preprocessing issues; Section 3.2 introduces the H- 



glCA model; Section 3.3 proves that H-glCA and GIFT coincide under noise-free settings; 
and Section |3.4| provides a measure of the functional homotopy for ICA-based networks. 
Section |4] provides a simulation study to demonstrate the effectiveness of H-gICA under 
homotopic, non-homotopic and mixed settings. This allows a comparison versus GIFT, 
demonstrating that by using intrinsic functional homotopy of the brain, H-gICA improves 
the power of locating the underlying brain networks. Section [5] provides the application 
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of the H-gICA on the ADHD-200 dataset, an open source dataset with 776 children and 
adolescents, while Section [6] contains a summary discussion. 

2 Background 

The observed data are denoted hj x = (xi, X2, - ■ ■ ,Xm)'^, an m-dimensional random vec- 
tor, and the true underlying signals by s = (si, S2, ■ ■ ■ , Sn)'^ , an n-dimensional transform 
of X. The problem is to determine a matrix W so that 

s = Wx. (1) 

A distinguishing feature of ICA components is that the elements of s are assumed to 
be independent in a generative factor analytic model, instead of focusing on data-level 
uncorrelatedness. The ICA problem can be formulated using the following generative 
model for the data: 

X = As, (2) 

where x is the observed m-dimensional vector, s is the ra-dimensional (latent) random 
vector whose components are assumed mutually independent, and A is a constant m x n 
matrix to be estimated. If it is further assumed that the dimensions of x and s are equal, 
i.e., m = n, the estimate of W in ([T| is then obtained as the inverse of the estimate for 
matrix A. 

Because of the assumption that A is invertible, a first stage dimension reduction is 
required. Thus the model can be more accurately stated as Kx = As, where 1^ is a 
dimension reduction matrix of size n x m. Normally, this is reached via a singular value 
decomposition of the demeaned data. Henceforth, we will assume that the data vectors 
are n dimensional. Also note that the data are typically de-meaned prior to analysis. 
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2.1 FastICA 



The previously mentioned dimension reduced data input into an ICA algorithm are mean 
zero and uncorrelated, and thus Gaussian distributional assumptions provide little further 
insight to linear reorganizations. This, along with the idea that linearly mixed iid noise 
tends to appear Gaussian, has lead researchers to design algorithms that search for an 
optimized W that maximizes non-Gaussianity. Such non-Gaussianity of the independent 
components is necessary for the identifiability of the model shown in Equation ^ ( Comon 



1994). Two common measures of non-Gaussianity used in ICA are kurtosis (Oja et al. 



2001) and negentropy (Pillai et al., 2002) 



Hyvarinen (1999) proposed a fast fixed-point algorithm known as FastICA. Following 



Hyvarinen's treatment of ICA, mutual information (Oja et al. , 2001), a natural measure 



of the dependence between random variables, can be written as 



when the variables i = 1, . . . ,n, are uncorrelated. Here y = {yi, ■ ■ ■ ,yn)'^ and J is 



negentropy, a measure of non-Gaussianity (Comon, 1994). Thus, W is determined so that 



the mutual information of the independent components, Sj, is minimized. As mentioned 



by Hyvarinen (1999), this is approximately equivalent to finding directions in which the 



negentropy is maximized. 

For any random variable, yi, with zero mean and unit variance, an approximation 
to the negentropy is proportional to [E{G{yi)} — E{G{Z)}]'^ for appropriately chosen 
G where Z is a standard normal. As an example, if G{y) = y^, one results in ICA 
methods that approximate minimum kurtosis algorithms. With this approximation, it 
can be shown that the problem is equivalent to finding argmax^. Y17=i Jci'^i) under the 
constraint that E{{w^x){w^x)} is with j ^ k and 1 when j = k. Here Jg is given 



by Jg{w) := [E{G{w'^x)} - ^{^(Z)}]^. The FastICA algorithm that we employ uses 
fixed-point iterations for maximization. 



2.2 Group ICA 

The extension of ICA to group inferences provides common independent components 
across subjects, which allows identification of putative common brain networks for the 
whole group. As in most fMRI studies, spatial independence is assumed in our group 
ICA model, since it is well-suited to the sparse distributed nature of the spatial pattern 



for most cognitive activation paradigms (McKeown et al. , 1997). However, instead of the 



entire brain, the spatial independence assumption in H-gICA is made for voxels within a 
single hemisphere. To use the information of brain functional homotopy, the fMRI data 
are registered to symmetric templates and thus the number of voxels are equal in the 



left and right hemispheres. Similarly to in Calhoun et al. (2001b), a dimension reduction 



using PCA is applied to each hemisphere of each subject. Another PCA step is then 



performed on the concatenated data matrix, which is not technically necessary (Eloyan 



et al. , 2013) though is the current standard practice for group ICA. 



3 Methods 

3.1 Preprocessing 

All group ICA methods for rs-fMRI require brains to be registered via a deformation al- 
gorithm to a canonical brain image, referred to as a template. This allows the assumption 



of spatially common brain networks elaborated on in Section 2.2 While our approach 
is applicable to any standard preprocessing, we do assume that the template is bilater- 
ally symmetric along the mid-sagittal plane, contrary to all human neuroanatomy. The 
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creation and improvement of such templates is beyond the scope of this work. However 
briefly, a simple strategy would require an accurate hemispheric model, a deformation 
to force the mid-sagittal plane to be perfectly planar, them flipping one hemisphere and 
pasting it together with itself to form the template. Fortunately, symmetric template 
brains exists; we employ one create by the International Consortium for Brain Mapping 
(ICBM) ( [Fonov et al. , 2011b, 2009a). An easy solution for preprocessing is to take already 
processed data, registered to a non-symmetric template, and apply a second registration 
to a symmetric template. As such, the method can be applied to existing processed data 
easily. Without loss of generality, it is assumed that the mid-sagittal plane is parallel to 
one of the image dimensions. 



3.2 Homotopic Group ICA 

nomotopic group ICA is the simple idea of having each subject contribute two fMRI 
images, one for their left hemisphere and one for their right, to the gICA approach. The 
benefit of the approach is the reduction in data level noise by (effectively) doubling the 
number of subjects and halving the number of voxels. Furthermore, it directly utilizes 
the information that many brain networks are apparently bilaterally symmetric. 

Suppose there are N subjects indexed by i = 1, 2, . . . , with an fMRI time series 
comprised of T scans. Further, suppose that each subject's fMRI image contains 2V 
voxels with V in each hemisphere for each time point. Separate the two hemispheres 
with two matrices of size T x V, labeled Fi^i and -Fj,2- Assume these have column sums 
zero (they are demeaned over time) and that one side has been appropriately flipped 
to geometrically match the other. Conceptually, let Fi{a,b,c,t) be a subject's image 
represented as a 4D array. Here, a indexes the left to right dimension where a = 1, . . . , 2K, 
b indexes the anterior/posterior dimension, c indexes the superior /inferior dimension and 
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t indexes scan. Suppose that the mid-sagittal plane is defined as a — K. Then, the 
4D array defined as Fi{a,b,c,t) for a < K defines one hemisphere's image series, while 
Fi{2K — a, 6, c, t) defines the other, oriented in a geometrically similar manner. The left, 
Fi^i, and right, -Fi,2, hemispheric data matrices are obtained from the associated left and 
right hemispheric images, having vectorized the spatial dimension. 

Assume a principal component dimension reduction has been applied and, abusing 
notation, that Xij of size Q x V denotes the dimension reduced data matrix of the left 
(j = 1) and right {j — 2) hemispheres of subject i. The column v of Xij continues to 
represent voxel v in hemisphere j. The rows of Xjj are the PCs, which are indexed by 
i = l,2,...,Q. 

Xij{t, v) represents row t, column v of X^j with the same convention applied to other 

vectors and matrices. Assuming that a group ICA decomposition implies: Xij{t,v) — 
Q 

^ Aij{t,q)S{q,v), for alH = 1,2, . . . ,iV and j = 1,2. This in turn implies that the 

9=1 

spatio-temporal process Xij{t,v) can be decomposed to a hemisphere specific time se- 
ries Aij{t,q) and a subject-independent spatial maps, S{q,v). This is equivalent to the 
following group ICA model: 

X = MS. (3) 

Here X — [X^, X^]^ is the 2NQ x V group data matrix from left and right hemi- 
spheres, where Xj := j = 1,2, which formed by concatenating 
N subjects' data in the temporal domain. S is Q x V matrix containing Q statistically 
independent spatial maps in its rows. M = [Aff , M'^]'^ is the 2NQ x Q group mix- 
ing matrix, where Mj — [A^^-, A^^-, • • • , Ajfj]"^ is the NQ x Q submatrix corresponding 
to hemisphere j concatenating the mixing matrix of the A'^ subjects. In the context of 
fMRI, the S{q,-), q — 1,2, ... ,Q are spatial maps that are the putative brain networks 
to be estimated. The mixing matrices are assumed non singular, and thus we can define 
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In the H-gICA model, the independent components are assumed to be common across 
subjects and hemispheres, while how they mix to produce the signal can differ among 
both subjects and hemispheres. Of course, the true mixing matrices are not observed and 
the FastICA algorithm is used to obtain the estimates. 



3.3 Connection with gICA 



In this section H-gIC A is linked to the most commonly used group ICA approach ( Calhoun 



et al. , 2001b). Let X = [Xi, X2] be the NQ x 2V dimension reduced group data matrix, 
where Xj = [Xjj, X^ j, ■ ■ ■ is same as in Section 3.2, Here, we are simply 



pasting the hemispheric data back together. The standard gICA approach decomposes 
X as X = MS. Here, S is a Q x 2V matrix containing Q spatial maps in its rows. 
And M = [Al , A2 , ■ ■ ■ , Atv ]^ is the NQ x Q group mixing matrix, where Ai is the 
Q X Q submatrix corresponding to subject i. Under the assumption that Ai, are of full 
rank, defines Wi = A^ ^, which can also be estimated by the FastICA algorithm. The 
following theorem, shows that if the actual sources are truly homotopic and noise free, 
H-gICA and gICA provide the same result when using the estimation method of FastICA. 
The proof of the theorem is given to the Appendix. 

Theorem 3.1. Suppose the actual sources are truly homotopic and noise free. The number 
of estimated ICs is Q <^V. Denote the FastICA estimate of S to be S and the estimate 
of S to be S. Then for j G {1, 2}, we have S = [S, S]. 



Theorem |3.1| shows that when there is no noise, H-gICA and GIFT are the same for 
the homotopic signals. However, of course, noise exists in most real data sets. Section |4] 
deals with H-glCA's ability to improve locating underlying sources when noise is added 
to the data. 
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3.4 Measures of Functional Homotopy 

In rs-fMRI, ICA output can be used to obtain measures of network connectivity, which 



Gao et al. 


2012 


)■ 


Joel et al. 


(2011) 



point out that the functional connectivity between any two brain regions may be due 
to within network connectivity (WNC) and between network connectivity (BNC). They 
emphasize the importance of interpreting such connectivity, ostensibly measured by the 
correlation and variance of the temporal mixing matrices. 



Similarly to Joel et al. (2011 ) the following defines a measure of subject- and network- 
specific functional homotopy for the fc*^ ICA based network of subject i: 

if.(A:) = Cor(AgUg), (4) 

where Aij {i = 1,2, . . . , N, j = 1, 2) is the mixing matrix of the i^^ subject corresponding 
to hemisphere j and A-^-^ is a vector of the time course modulating spatial map k. Note, 
typically this requires a back transformation from PCA space. Here Hi{k) measures 
the spontaneous activity of the /c*'^ network between left and right hemispheres. The 
estimation of Hi{k) in H-gICA is given by replacing Ai^i and Ai2 with their estimated 
values in (|4]). Similarly the group level functional homotopy for the k^^ ICA based network 
can be define as: 

H{k) = CoT{M[^\Mi^^), (5) 

where Mj is defined in (|3| as the submatrix corresponding to hemisphere j concatenating 
the mixing matrix of the subjects and A4"j*^^ is the t*'^ element of the time course 
modulating spatial map k. 
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4 Simulation Results 



4.1 A Simple Example 



The following simple example is given to initially demonstrate the effectiveness of H- 
glCA. Suppose the number of subject is = 3 and the number of underlying sources is 
Q = 3. For each subject i, the T x 2V data matrix = ^4,2] is from the model 

Xj . = AjS with T = 3 and 2V = 10^. We further assume that 



^ -1 -5 2 ^ 



-3 2 



y 5 3 -5 y 



^ -1 -1 2 ^ 



^-4 5 y 



^-3 3 -5^ 



-1 -1 



-2 -1 



(The elements in Ai, A2 and A3 are randomly picked from Unif{— 10, 10).) 

As shown in Figure [l| three different source matrices (S) are generated: one is sym- 
metric; another is only in one hemisphere; and the third has two asymmetric blocks of 
activated voxels. These results are shown in the bottom row of Figure [T| highlighted by 
H-glCA's separation of the three source matrices and its meaningful prediction for all of 
them. 



4.2 2D Simulation 

In order to illustrate the performance of the proposed method, simulation studies with 
10, 000 voxels in 2D spaces were conducted. Both the homotopic and non-homotopic 
settings are used in the study. The results are compared with the commonly used group 
ICA algorithm without consideration of the left and right hemispheres. 

Similar to above, the number of subjects is A^ = 3 and the number of underlying 
sources is Q = 3. The data are generated by the ICA model Xi^. = AiS with T = 3 and 
2V = 100^. Ai has the same value as in the toy example for 2 = 1, 2, 3 is also assumed. 
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□ B □ 

r ri 

Figure 1: Result of the first example. The top row shows the true sources in three 
different types: perfectly symmetric; only present in one hemisphere; differing in the two 
hemispheres. The small blocks in these plots are activated voxels, where the values come 
from uniform distribution. The first two elements of the second row are the true sources 
in the left hemisphere and the last element is left minus right. The bottom row consists 
of the independent components generated by H-glCA. 
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4.2.1 Case I: Perfect Homotopy 

Assume all true sources are perfectly homotopic. For each source, two blocks of voxels are 
symmetrically activated. In each loop of the simulation, values of these activated voxels 
are assigned by Gamma distributions. The heatplots of the sources are shown in Figure 
[2j The data are generated by the three sources in Figure |2] and the mixing matrices 



Source 1 Source 2 Source 3 




shape=7. rate=5 shape=8, rate=4 shape=9. rate=3 



Figure 2: True Sources in the Perfect Homotopy. The three sources are generated by 
Gamma distributions with different parameters. 

Ai, A2, and A3. Gaussian noise is then added to the data. The estimated components 
are standardized and subsequently compared with the corresponding standardized true 
sources. Next, the mean difference at each voxel is calculated. The simulation contained 
300 iterations in each run. Homotopic gICA results are compared with the commonly 
used group ICA algorithm in Table [1} where the noise is set to be mean zero and with a 
standard deviation equal to 5. 

As seen in Table [T| in the setting of symmetric sources, the mean errors of H-gICA 
are smaller than that of ordinary glCA for all the three sources. Figure |3] shows the 
mean of the voxel mean difference with different settings for the noise. Again, H-gICA 
works consistently better when noise exists and is the same as ordinary glCA in noise-free 
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Table 1: The two rows compare the mean of the voxel mean difference of H-gICA and 
ordinary gICA in the symmetric setting. 





Source 1 


Source 2 


Source 3 


H-gICA 
gICA 


0.414 (~ 0.001) 
0.507 (~ 0.001) 


0.412 (~ 0.001) 
0.496 (~ 0.001) 


0.423 (~ 0.002) 
0.433 (~ 0.001) 



settings, which was been proven in Theorem 3.1 



Figure |4] compares the ICs estimated by H-gICA and the ordinary glCA. As we can 
see, comparing with the ordinary gICA, H-gICA has a consistently better estimation of 
the ICs when noise exists. 

4.2.2 Case II: Perfect Lateralization 

True sources are now presumed to be only present in one hemisphere (perfect lateralization 
of the brain network), and without loss of generality, the left hemisphere is selected. 
Similar to Case I, a Gamma distribution was used to generate the value of the activated 
voxels and, again, 300 iterations were run. Figure [5] shows the heatplots of the three 
true sources in the first iteration. Gassian noise was added to the data generated by 
the three sources and the results are compared with ordinary glCA. Note that, since the 
independent components provided by H-gICA contain only one hemisphere by design. 
Thus, the sole applicable comparison is the true sources in the left hemisphere. Table |2] 
gives these results, and show that the mean error of H-gICA is less than that of gICA for 
Source 2 while larger than the mean error of gICA for Source 1 and Source 3. In summary, 
H-gICA is competitive with normal g-ICA when the sources are lateralized with the caveat 
that H-gICA does not provide lateralization information on which hemisphere the network 
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standard error of noises 



Figure 3: The mean of the voxel mean difference increases with the standard deviation 
of the noise. The two methods are the same under noise-free settings and when the noise 
exists, H-gICA works better than ordinary glCA. 
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Figure 4: The ICs estimated by H-gICA and glCA. The first column shows the true 
sources. The other columns, from left to right, show the estimated ICs when the noise 
increases. The 1^*, 3'"'^, and 5*^ rows are for H-gICA and the 2nd, 4*^* and 6*^ are for glCA. 




Source 3 




25 50 75 100 

shape=7, rate=5 shape=8, rate=4 shape=9, rate=3 

Figure 5: True Sources in Non-homotopic Setting. All three true sources are only in one 
hemisphere. The values of the activated voxels are generated by Gamma distribution with 
different parameters. 
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resides in. We do note, however, that this information is contained in the temporal mixing 
matrix, just not in a form easily displayed as an image. 

Table 2: Comparison of the H-gICA results versus ordinary gICA results when the true 
sources are only in one hemisphere. The table provides the mean and standard deviation 
of the voxel mean difference with the true sources for the 300 iterations. 





Source 1 


Source 2 


Source 3 


H-gICA 
gICA 


0.508 (~ 0.002) 
0.502 (~ 0.003) 


0.494 (~ 0.001) 
0.533 (~ 0.002) 


0.431 (~ 0.002) 
0.312 (~ 0.003) 



4.2.3 Case III: Mixture of Lateralized and Homotopic Networks 

For this setting, a mixture of two different types of sources is introduced. As shown 
in Figure [6} the first two sources are homotopic while the third source covers different 
regions in the two hemispheres. The estimated ICs for the two homotopic sources are in 
Figure [7} which illustrates that the effectiveness of estimating the homotopic sources is 
not impacted by adding non-homotopic sources. 



4.3 Flag Example 

In this example, the actual sources will be the gray-scaled flags of the USA, Canada, the 
European Union, China and Russia. As shown in Figure [8| three of them are symmetric 
(Canada, the European Union and Russia) and two are not (USA and China). Similar 



as in Section 4.2, the data are generated by the ICA model with a fixed mixing matrix, 
Gaussian noise was then added. The results are shown in Figure^ where the I''*, 3'''^ and 



5*^* rows are the estimated symmetric sources extracted by H-gICA and the 2ndi 4*'^ and 
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Source 3 




25 50 75 100 

shape=9, rate=3 shape=9. rate=3 shape=7, rate=5 

Figure 6: True Sources in Mixed Setting. The first two sources are homotopic while the 
third one varies in the two hemispheres. 
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Figure 7: The homotopic ICs estimated by H-gICA and ordinary glCA. The first column 
shows the true sources. The other columns, from left to right, contain the estimated ICs 
with increasing noise. The l*** and 3'"'^ rows are results of H-gICA and the 2nd and 4*'* are 
of ordinary glCA. 
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6*'* rows are the estimated symmetric sources extracted by ordinary glCA. As we can see, 
for all of the three symmetric sources, H-gICA provides clearer estimation as the noises 
increased. Moreover, leakage that is apparent in the g-ICA is not apparent in H-glCA. 




Figure 8: Actual sources of the flag example. The original flags images were taken from 
Wikipedia. 



5 Application to the ADHD-200 Dataset 

Application of the H-gICA method is illustrated using the ADHD-200 dataset, one of the 
largest and freely available resting state fMRI datasets. The data combines 776 resting- 
state fMRI and anatomical scans aggregated across eight independent imaging sites, 491 
of which were obtained from normally developing individuals and 285 from children and 
adolescents with ADHD (ages: 7-21 years old). We view this analysis as largely a proof 
of principle in applying the method and defer thorough investigations of ADHD for later 
work. 

This particular analysis focused on 20 subjects picked from the ADHD-200 data set. 
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Figure 9: The estimated homotopic sources of H-gICA and ordinary glCA. The five 
columns, from left to right, contain the estimated ICs with increasing noise. The 1**, 3'^'^ 
and 5*^ rows shows the estimated sources extracted by H-glCA while the 2nd, 4*^* and 6^^ 
rows shows the estimated sources extracted by ordinary glCA. 
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Data were processed via the NITRC 1,000 Functional Connectome processing scripts 



(Mennes et al. , 2012). In summary, images were slice-time corrected, deobliqued, skull 
stripped, smoothed and registered to a 3 mm3 MNI template. The data were then regis- 
tered to ICBM 2009a nonlinear symmetric templates generated by the McConnell Brain 



Imaging Centre (Fonov et al. , 2009b, 2011a). Each fMRI scan contains 99 x 117 x 95 



1, 100,385 voxels measured at 176 time points. Figure [10] and Figure [TT] are the QQ plot 
and scatter plots of the estimated sources extracted by ordinary gICA in the left and right 
hemispheres. Most of the estimated sources are close to the 45° line, which suggests that 
the marginal distributions of left and right hemispheres are similar. H-gICA can benefit 
from this apparent homotopy. 




Sources in Left Hemisphere 



Figure 10: QQ plot of the 15 sources extracted by glCA. 



Our procedure is as follows. Each fMRI scan was separated into left and right hemi- 
spheres. Thus, each hemisphere contained 49 x 117 x 95 = 544,635 voxels. Similar to 



standard group ICA (Calhoun et al. 2001b), a dimension reduction using PCA was ap 



plied to each hemisphere of each subject. 15 PCs are obtained for each hemisphere. A 
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Figure 11: The scatter plot of the 15 sources. 



group data matrix was generated by concatenating the reduced data of both hemispheres 
of the 20 subjects in the temporal domain. Thus, the aggregated matrix has dimension 
2NT X V, where N = 20, T = 15, and V = 544, 635. Our algorithm of homotopic 
group ICA is then applied on this matrix. Fifteen estimated independent components 



are postulated by H-glCA. As shown in Figure 12, out of the 15 components, several 



brain networks were found including: the visual network 12(a) , the default mode network 



12 c 



the auditory network 12(e) , and the motor network 12(g) Compared with the ICs 
obtained from ordinary gICA, shown in |12(b) 12(d), 12(f) and 12(h), H-gICA improves 
the estimation of all of these sources by yielding substantially more clearly delineated 
networks. 

The approach of H-gICA also allows us to calculate the brain functional homotopy of 
each brain network. To compare the brain functional homotopy of ADHD and typical 
developed children, we choose 20 ADHD subjects and 20 subject-macthed controls. The 
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subjects and controls were matched in gender and age. Via Equations ^ and ([s]), the es- 
timated functional homotopy of four networks (visual, default mode, auditory and motor) 



are shown in Figure [T3| As we can see, the functional homotopy of ADHD children tends 
to be lower in both visual networks and the auditory network. These represent meaningful 
leads on the exploration of homotopic network relationships and disease, though we leave 
a full exploration to later work. 



6 Discussion 

In this paper we present a new group ICA method called homotopic group ICA (H-gICA). 
Similar with ordinary group ICA methods, H-gICA can analyze data for multiple sub- 
jects and estimate common underlying IC's across individuals. By concatenating the fMRI 
data of the two hemispheres, H-gICA effectively doubles the sample size. It improves the 
power of finding the underlying brain networks by rearranging the data structure and 
utilizing the known information of synchrony between bi-laterally symmetrically oppos- 
ing inter-hemispheric regions. Both the simulation study and the application on ADHD 
200 data show that H-gICA is preferable to ordinary gICA when the data are homo- 
topic, increasingly so as noise increases. Moreover, H-gICA remains preferable to gICA at 
estimating homotopic sources, even in the presence of non-homotopic sources. Effective- 
ness was demonstrated by application on the ADHD-200 dataset. Several brain networks 
were found and clearly represented in smoother, more clearly delineated, contiguous vol- 
umes than ordinary glCA. The main networks found included visual networks, the default 
mode network, the auditory network, as well as others. Moreover H-gICA enables certain 
measurement of the functional homotopy of the underlying functional networks. This po- 
tentially offers the opportunity to analyze the relation of the brain functional homotopy 
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• Normal 

• ADHD 




Figure 13: Comparison of functional homotopy of ADHD and normal developed children. 
Each column represents a network (visual, default mode, auditory and motor). Each pair 
(subjects and controls) are connected via a grey line. The left end points of the grey 
lines measure the functional homotopy of the control and the right end points are for the 
ADHD subjects. The black lines represent the group level functional homotopy for the 
four networks. 
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between the left and right hemisphere of the brain with diseases. 
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Appendix 



Proof of Theorem 13.1 



Proof. Since the true sources are truly homotopic and noise free, we have: 

where Xjj, j = 1,2, are the data of left and right hemispheres after dimension reduction 
by PCA. This is equivalent to 

Xi = X2- 

Without loss of generality, we assume Xj {j = 1, 2) are demeaned i.e. the row means of 
Xj are all 0. Assume the singular value decomposition of the matrix Xi/\/V is 



Xi/W = u^v 



T 

5 



where U is of dimension NQ x {NQ — 1) which consists of the left singular vectors of 
Xi, S is a diagonal matrix of dimension {NQ — 1) x {NQ — 1) with the singular values, 
and V is of dimension V x {NQ — 1) which consists of the right singular vectors. Both 
U and V are orthogonal. Note that, since X is not full rank, only NQ — 1 singular value 
are non-zero and thus only NQ — 1 singular vectors are estimated. Then we have 

x/W = [{xjWf, {x^/Wff 
= [{uj:v^f, {UEV^'ff 

and 

X/W =[{Xi/Vv),{X2/Vv)] 

= [(C/SV^), (C/SV^)] 

= uj:[v^,v^]. 
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Thus we can define 




(6) 



to be the pre-whitening matrix that projects data onto the principal components: 

Z = KX and Z = KX, 

where Z = \fVV^ and Z = {\fVV^ , \lVV^\ are the whiten data for H-gICA and GIFT 
respectively. Clearly, if we assume the random variables z and z are taking values from 
the columns of Z and Z^ then we will have: 



Suppose we have the same initial value for the vectors lOp for all p e {1, 2, • • • , Q}. 
No matter how the contrast function, G, is defined for each fixed w we have: 



E{zz^\ = I and E[zz'^] = I. 



E{zgiw^z)} = i^Z(,^Mw,^Z(,^;)) 



v=l 



V 



v=l 




-^Y.~^{:v)g{w^~Z{-,v)) 



v=l v=l 



v=l 



= E{zg{w'^z)} 



where g is defined as the derivative of the contrast function G. Similarly, 



E{g'{w^z)}w ^ E{g'{w'^z)}w, 



where g is the derivative of G. So we have: 



E{zg{w'^z)} - E{g'{w'^z)}w = E{zg{w'^z)} - E{g'{w'^z)}w 
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Thus, following the FastICA algorithm (Hyvarinen, 1999), it is easy to see that the esti 



mates of Wp, p = 1, 2, ■ ■ ■ , m, will be same for these two approaches. By the definition of 
K and ^ in (|6|, we will have: 



and 



1 



s = wx = -[w,w][xi,x'2Y = wxi 



S = WX = W[Xi,X2] = W[Xi,Xi] 
Thus follows the result of Theorem 13.11 



□ 
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